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Abstract 

We present results of a series of Monte Carlo simulations investigating the 
imprint of a central intermcdiatc-mass black hole (IMBH) on the structure of a 
globular cluster. We investigate the three-dimensional and projected density pro- 
files, and stellar disruption rates for idealized as well as realistic cluster models, 
taking into account a stellar mass spectrum and stellar evolution, and allowing for 
a larger, more realistic, number of stars than was previously possible with direct 
N-body methods. We compare our results to other N-body and Fokker-Planck 
simulations published previously. We find, in general, very good agreement for 
the overall cluster structure and dynamical evolution between direct iV-body 
simulations and our Monte Carlo simulations. Significant differences exist in the 
number of stars that are tidally disrupted by the IMBH, which is most hkely an 
effect of the wandering motion of the IMBH, not included in the Monte Carlo 
scheme. These differences, however, are negligible for the final IMBH masses in 
realistic cluster models as the disruption rates are generally much lower than for 
single-mass clusters. As a test to observations we model the cluster NGC 5694, 
which is known to possess a central surface brightness cusp consistent with the 
presence of an IMBH. We find that not only the inner slope but also the outer 
part of the surface brightness profile agrees with corresponding observations. 
This provides further evidence for the possible existence of a central IMBH in 
NGC5694. 

1. Introduction 

As recently as 10 years ago, it was generally believed that black holes (BHs) occur 
in two broad mass ranges: stellar {Mbh — 3 — 2OM0) , produced by the core collapse of 
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massive stars, and supermassive {Mbh ~ 10^ — lO^'^M©), believed to have formed in the 
centers of galaxies at high redshift and grown in mass as the result of galaxy mergers 
(see, e.g., Volonteri, Haardt & Madau 2003). However, the existence of BHs with masses 
intermediate between these two ranges could not be established by observations up until 
recently, although intermediate mass BHs (IMBHs) were discussed by theorists more than 
30 years ago; (see, e.g., Wyller 1970). Indirect evidence for IMBHs has accumulated over 
time from observations of so-called ultraluminous X-ray sources (ULXs), objects with 
fluxes that exceed the angle- averaged flux of a stellar mass BH accreting at the Eddington 
limit. An interesting result from observations of ULXs is that many, if not most, of them 
are associated with star clusters. It has long been speculated (e.g., Frank & Rees 1976) 
that the centers of globular clusters (GCs) may harbor BHs with masses ~ IO^Mq. If so, 
these BHs affect the distribution function of the stars, producing velocity and density cusps 
(Bahcall & Wolf 1976). While the detection of ULXs can only give indirect evidence of 
the presence of IMBHs, observations of cuspy velocity profiles would make it possible to 
directly determine the BH mass. However, the radius of influence of an IMBH, defined as 
the radius where the orbital velocity around the BH equals the velocity dispersion of the 
cluster, is very small. For example, at a distance of lOkpc, a lO^M© BH would influence 
orbits within 1" , making observations very challenging. 

Studies of the surface density proflle of GCs offer a complementary method of 
constraining the effects of an IMBH on the host GC stars. A recent study by Noyola & 
Gebhardt (2006) obtained central surface brightness proflles for 38 Galactic GCs from HST 
WFPC2 images. They showed that half of the GCs in their sample have slopes for the 
inner surface brightness proflles that are inconsistent with simple isothermal cores, which 
may be indicative of an IMBH. However, it is challenging to explain the full range of slopes 
with current models. While analytical models can only explain the steepest slopes in their 
sample, recent N-body models of GCs containing IMBHs (Baumgardt et al. 2005), might 
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explain some of the intermediate surface brightness slopes. 

However, the disadvantage of current N-body simulations is that for realistic cluster 
models, which take into account stellar evolution and a realistic mass spectrum, the number 
of stars is restricted to typically less than ~ 10^ as these simulations require a large amount 
of computing time. However, many GCs are known to be very massive, with masses 
reaching up to 2 x 10^ resulting in a much larger number of stars one has to deal 
with when modeling these objects. In previous N-body simulations, such large-N clusters 
have been scaled down to low-N systems. Scaling down can be achieved in two ways (e.g. 
Baumgardt et al. 2005): either the mass of the central IMBH Mbh is kept constant and 
is decreased, effectively decreasing the total cluster mass Mc, or the ratio Mbh/Mq is 
kept constant, while lowering both Mbh and Mc- As both Mbh/Mc and the ratio of Mbh 
to stellar mass are important parameters that influence the structure and dynamics of a 
cluster, but cannot be held constant simultaneously when lowering A^, it is clear that only 
with the real A^ a fully self-consistent simulation can be achieved. Scaling becomes even 
more difficult once other physical processes are included into the simulations such as stellar 
evolution or stellar collisions. Using the correct number of stars in a dynamical simulation 
ensures that the relative rates of different dynamical processes (which all scale differently 
with A^) are correct. 



It is clear that in order to study the evolution of old globular clusters that might 
harbor IMBHs, more approximatee methods have to be employed. They fall roughly 



(Amaro-Seoane e 


t al. 


2004; 


Takahashi 1997: 


Giersz & Spurzem 


1994: M 


urphv et al. 


1991; 


Cohn & Kulsrud 


1978 


Liehtman & Shapiro 


1977; 


Bahcall & WoU 


1977 


1, and Monte Carlo 


(MC) methods t 


lat use a 


particle based approach (see, e.g., 


Fret 


reau & Rasio (200 


7); 


Freitae & Benz ( 


2002. 


2001 


): Shapiro & Marchan 


; (1978): Henon 


(1971 


) and references 



- 5 - 



therein). While the former methods ha ve been successfully used to study the effect of 



a massive BH on the cluster structure (Amaro-Seoane et al 



Cohn &: Kulsrud 



1978 



Lightman fc Shapiro 



2004 



Murphy et al 



1991 



19771 ). the model clusters were highly idealized, 



consisting only of equal-mass stars, and did not incorporate stellar evolution. Although 
including additional physical processes is not impossible, it remains nevertheless highly 
non-trivial for these methods. The MC method, on the other hand, relies on a star-by-star 
description of the cluster and has, therefore, the great advantage that additional processes 
are easily incorporated. 

Our group has been developing over many years a state-of-the-art MC code, which 
treats many relevant processes in sufficient detail, making direct comparison with GC 



observations feasible ( jFregeau fc Rasio 



120071 and references therein). This paper is the sixth 
paper in a series studying the fundamental aspects of cluster dynamics using this code. 
Here we will describe the changes we made to our code in order to incorporate the effect of 
a massive central IMBH and carry out comparison runs with idealized models as well as 
more realistic cluster simulations published previously in the literature. In ^ we briefly 
describe our method and the changes we made to the MC code. We validate our code 
by comparing our results to previously published results in the literature, using idealized 
cluster models (§!]), as well as more realistic ones that include stellar evolution ([|5]). In 



§5.21 we present 
observations of 



surface brightness 



j roflles from our large- runs and compare them to 



Noyola &: Gebhardt 



(120061 ). We conclude in 



2. Previous Work of Globular Cluster Evolution with IMBHs 

The dynamical effect of an IMBH on the surrounding stellar system was flrst described 
by Peebles (1972), who argued that the bound stars in the cusp around the BH must 
obey a shallow power-law density distribution to account for stellar consumption near the 
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cluster cente r. Analyzing the 
distribution 



bkke r-Planck equation in energy space for an isotropic stellar 



Bahcall fc Wolj ( 19761 ) obtained a density profile with n{r) oc r^'^/'^, which is 
now commonly referred to as the Bahcall- Wolf cusp. The extension of the cusp solution is 
given by the radius of influence of the BH, rj, which is defined as 

GMbh 



1 



where G is the gravitational constan t and a the velocity dispersion of the core. As shown 



later by lShapiro fc Lightmanl (jl976l ). such a solution can be readily obtained using simple 
scaling arguments. The key is to realize that in the region delimited by the tidal radius, r^, 
within which stars get tidally disrupted, and within which the stars are bound to the 
black hole, the net energy diffusion timescale, is proportional to the local relaxation 
time, tr, which is the shortest timescale on which any physical quantity can be transported. 
Furthermore, the quasi-steady-state of the cusp region is characterized by a dynamic 
equilibrium, with a constant net energy flow into the core region that should scale as 
n{r) r^E{r)/tu, where E{r) is the mean specific energy at radius r, so E{r) ~ GMbh't'^ 
in the cusp region. Setting then tu ty. ^ a^/^/ [G'^m?n{r)\ , a ~ ^/GMBIpr~^, and simple 
substitutions immediately lead to n(r) oc r 

The formation of such a cusp has been confirmed subsequently by numerical 



studies. M ost of them emp 



it directly flBahcall fc WoU 



oyed methods based on the Fok 



1976 



Lightman fc Shapiro 



^ er-Planck equation, so 



1977 



Cohn fc Kulsrud 



or inc 


irectly, either using t 


1978 




Freitag & Benz 


2002) 


( Amaro-Seoane et al. 


2004) 



1978), 



vmg 



2002) or a fluid-dynamical approach based on velocity moments 
20041 ). Being derived from the Fokker-Planck equation, all these 
methods share essentially the same set of underlying assumptions: (i) the cluster potential 
has spherical symmetry; (ii) the cluster is in dynamical equilibrium at all times; (iii) the 
evolution is driven by diffusive 2-body relaxation. Through direct A^-body simulations that 
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do not rely 



solution of 



on any a prior i assum ptions, iBaumgardt et al.l (l2004a[ ) confirmed the cusp 



Bahcall fc Wolj (119761 ) . therefore also providing important justification to the 



Fokker- Planck approach. 



Based on the cusp solution, 



Frank fc Reed (119761 ) calculated the stellar disruption 



rate taking into account that stars inside a critical radius, Vcrit, are efficiently accreted by 
the black hole as they diffuse quickly into low angular momentum orbits with periastron 
distances, Vperi, smaller than rj. As for a given radial position, the velocity vectors that lead 
to orbits with rperi < rt, form a cone, these orbits are also called loss-cone orbits. Outside 
Vcrit, stars are always able to leave the loss-cone during one orbital period due to two-body 
relaxation while inside Vcrit the orbital changes are smaller so that stars on loss-cone orbits 
are more likely to reach the tidal radius and get disrupted before they have a chance to get 
scattered out. Consequently, inside Vcrit the loss-cone should b e nearly empty while outside 



Vcrit it always remains full. In addition. 



Frank &: Reed (119761 ) argue that the disruption 



rate is mainly given by the cluster conditions at Vcrit- This is because, on the one hand, 
for r > rcrit the fraction of stars populating the loss-cone decreases as the loss-cone angle 
decreases with increasing radius, while for r < r^it the total net flux of stars, and therefore 



the fl ux of stars that c an diffuse into the los s -cone , decreases rapidly (ILightman fc Shapiro 



19771 ). Calculations by 



Amaro-Seoane et al 



( I2OO4I ) conflrm the loss-cone picture, showing 



that, for a cluster in dynamical equilibrium, the disruption rate is strongly peaked at Vcru, 
and the fraction of stars on loss-cone orbits is ra pidly approaching zero fo r r < Vcru, while 
for r > Tcrit it is always close to one. Similarly, 



Baumgardt et al. 



( j2004al ) find generally 



good agreement between the disruption rates in their simula tions and the disruptions rates 



based on the approximate expression of 



Frank fc Reed (119761 ) (their equation 22) and using 



the cluster conditions at rcrit from their A^-body models. 



While earlier work on the dynamics of clusters with IMBHs mainly focused on the 
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equilibrium state of the cusp surrounded by a static isothermal core, IShapird (119771 ) 
considered the effect of an IMBH on the global evolution of the cluster. Using a homological 
model for the dynamical evolution, he calculated the core size in response to evaporation 
of high-velocity stars and tidal disruption of stars tightly bound to the IMBH. While stars 
that leave the cluster by evaporation carry away very little energydriving core contraction 



( ISpitzer fc Saslaw 



19661 ). which would ultimately lead to core collapse in the absence 
of an IMBH, the tidal disruptions close to the central IMBH that remove stars with 
highly negativ e specific energies, provide an energy source that causes the core to expand. 



Shapiro! ( 119771 ) shows that for low initial IMBH masses and large initial core radii, stellar 
evaporation first dominates and drives core contraction until, due to the increasing core 
density, the tidal disruption rate becomes large enough to reverse core collapse. The time 
of this reversal roughly coincides with the time of core collapse for the cluster without 
IMBH. Tidal disruptions then drive the re-expansion of the core, and the core size increases 
asymptotically to infinity. This expansion is a generic feature of a stellar system where 



region is very small compared to the cluster mass (iHenon 



Henon 


1965; 


Shapiro 


1977) 



The qualitative behavior of the c ore size evolution for lower mass IMBHs was 



later confirmee 



Freitag &: Benz 



by numerical studies (IMarchant &: Shapiro 



2002 



Amaro-Seoane et al. 



20041 ) 



1980 



Murphv et al 



Amaro-Seoane et al. 



1991 



(120041 ) in particular 



showed that t 



le core size increases asymptotically as oc t^^^, which was also predicted by 



Shapiro! (119771 ). This e xpansion also causes the d isruption rate to decrease with time as 



approximately oc t (Amaro-Seoane et al.ll2004l ). and will ultimately l ead to the com plete 



dissolution of the cluster as the outer stars are removed by tidal forces (jWielen 



197l|). For 



larger IMBH masses and small initial core sizes, tidal disru ptions will pre vent any initial 



core contracti on and the core e x pands from the beginning ( IShapiro 



calculated by 



Baumgardt et al. 



1977! ). This case was 



(l2004al ) using direct A^-body simulations, confirming that 
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core expansion starts almost immediately and follows a t^/^ power-law. 



Most of the studies mentioned so far considered the evolution of cluster containing 



a central massive black hole anc 
spectrum was fi r st con sidered by 



comprised of s t ars o: 



" equal mass. A stellar mass 



Bahcall fc WoliI (119771 ) extending their previous work in 



Bahcall fc Wolj (119761 ). They find that, due to mass- segregation, lower mass stars have 
shallower density profiles than more massive ones. For old globular clusters this means that 
the observable surface brightness cusp must be much shallower as the more massive dark 
stellar remnants are concentrated towards the center while the lower-mass main sequence 
stars that contribute most of the light are much less centrally concentrated. It follows 
that, although a cusp in the velocity and density profile provides strong evidence for the 
presence of an IMBH in a cluster, such cusps might not be easily detectable in real star 
clusters. Using d irect iV-body simulations including an initial stellar mass spectrum and 



stellar evolution, 



Baumgardt et al 



(l2004bl ) find, indeed, flat luminosity density profiles 



almost indistinguishable fr om a standard Ki n g pro file. Carrying out similar simulations 



but with Mbh/Mc < 1%, iBaumgardt et al.l (120051 ) find surface brightness cusps with 



power-law slopes ranging from a = —0.1 to a = —0.3. Based on these results they identified 
9 candidate clusters from the sample of galactic GCs of Noyola & Gebhardt (2006) that 
might contain IMBHs. 

Monte Carlo simulations of realistic clusters with central IMBH were mainly done 



in the context o: 



gala ctic nuclei. A recently developed and well tested code is that o f 



Henonl (Il97lh but 



Freitag Sz Bena (120021 ). Similar to our code, it is based on the method of 
is modified to evolve each star individually on a fixed fraction of the local relaxation time, 
as opposed to the original shared time-step scheme. In addition to the implementation of 
loss-cone physics, which we will describe in detail in ^ it also incorporates stellar collisions 
interpolating between results from detailed hydrodynamical simulations. Collisions between 
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sta rs is an impor t ant p hysical process in dense galactic nuclei. As already pointed out 



by iFrank &: ReesI (Il976l ). the radius Vcou outside which stellar encounters responsible for 
relaxation can be treated as elastic encounter can be larger than Vcrit for large Mbh and 
typical sizes and densities for galactic nuclei. Inside rcou stars cannot deflect each other 
significantly without colliding. As the relative velocity within Vcoii is larger then the escape 
velocity from the stellar surface the collision can be very disruptive and mig ht under certain 



conditions provide a significant source to fuel an active galactic nucleus (see 



2002 



Freitag fc Benz 



, and references therein). However, for conditions typical inside globular clusters stellar 
collisions are unlikely to play a significant role and are therefore not further considered for 
the present study. 

In addition to the formation of cuspy profiles, an IMBH influences the surface brightness 



profile by producing rather large cluster cores as measured b y the core-to- 



half 



ight radius, 
20071 ). The large 
le core 



such that larger cores are produced by more massive IMBHs (iHeggie et al. 
core sizes are simply a result of the energy flow from the central cusp regio n to t 
which causes the cluster to expand. Constructing generalized King rnodels 
including the effect of an IMBH and a stellar mass spectrum, iMiocchil (120071 ) finds that the 



(jKing 



19661) 



core size and the cusp slope are related such that clusters with larger slopes, s, have lower 
concentrations, c = log(rc/r(), where Vc is the core radius of the cluster. More specifically 
they find that s and c are related by 



11.6s - 4.85 < log ( 1 < -1.14c - 0.694 . 



where Mr is the cluster mass. Based on this criterion and data from 



(2) 



Noyola fc Gebhardt 



(120061 ) for s as well as the Harris catalog for c, they identified 7 candidate clusters that migh t 



contai n IMBHs with mass > 100 M©, with 4 of them also identified by 



(l2005h 



Baumgardt et al. 



In contrast to the rather shallow surface brightness cusps, the stronger cusp in the 
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stellar velocity dispersion appears to be a much better diagnostic to infer the presence of 
IMBHs in globular clusters. However, for globular clusters this signature turns out to be 
difficult to detect as any velocity dispersion measurement inside such a cusp has to rely on 
only a f ew bright stars for expe cted IMBH masses < 1000 Mq and typical globular cluster 



masses ( iBaumgardt et al. 



20051 ). These IMBH mass estimates are based on extrapolating 
the well known Mbh ~ <^ relation between central black hole mass and velocity dispersion 
in the central bulges of galaxies down to velocity dispersions typical for globular clusters 



lOkms ^). In that case Mbh sho u. 



collisional runaways by 



Giirkan et al. 



d be ~ 10'^ — 10^ Mq. Furthermore, simulations of 



(120041 ) show that the mass of the Spitzer unstable 



subcluster which provides the mass reservoir for forming the BH progenitor is ~ 10^ Mc, 
implying IMBH masses of at most a few 10^ Mq for typical cluster masses. However, it is 
important to point out that this does not mean that Mbh/Mc for an old globular cluster 
has to be always significantly less than 1%, as a cluster can later loose a substantial amount 
of mass due to tidal stripping in the Galactic field. 



Method and Initial Conditions 



3.1. Monte Carlo Method with IMBH 



Our MC code shares some important properties with direct A^- 



body methods, which is 



Freitag &: Benz 



200l|) 



why it is also regarded as a randomized A^-body scheme (see, e.g. 
Just as in direct A^-body codes, it relies on a star-by-star description of the GC, which 
makes it particularly straightforward to include additional physical processes such as stellar 
evolution. Contrary to direct A^-body methods, however, the stellar orbits are resolved on 
a relaxation time scale tr, which is much larger than the crossing time tcr, the time scale on 
which direct A^-body methods resolve those orbi ts. The spec i fic irn plementation we use for 



our study is the MC code initially developed by 



Joshi et al. 



(I2OOOI ) and further enhanced 



and improved by iFregeau et al.l (120031 ) and lFregeau &: Rasiol (120071 ). The code is based on 
Hnon's algorithm for solving the Fokker-Planck equation. It incorporates treatments of 
mass spectra, stellar evolution, primordial binaries, and the influence of a galactic tidal 
field. 

The e ffect of an IMBH on t he stellar distribution is implemented in a manner similar 



to that of iFreitag fc Bend (|2002| ). In this method the IMBH is treated as a fixed, central 
point mass while stars are tidally disrupted and accreted onto the IMBH whenever their 
periastron distances lie within the tidal radius, r^, of the IMBH, which is given by 



M, 



BH 



1/3 



(3) 



where i?* and M^, are the stellar radius and mass respectively. Stars are removed from the 
system and their masses are added to the BH as soon as their velocity vector, v, enter the 
loss-cone, 6lc^ approximately given by 

.GMBHTt 



However, as the star's removal happens on an orbital timescale one would need to use 
timesteps as short as the orbital period of the star in order to treat the loss-cone effects 
in the most accurate fashion. This would, however, slow down the whole calculation 
considerably. Instead, during one MC timestep a star's orbital evolution is followed 
by simulating the random-walk of its velocity vector, which approximates the effect of 
relaxation on the much shorter orbital timescale. The random-walk procedure goes as 
follows: 



1. After a gravitational encounter between two stars is calculated in the standard Monte 
Carlo fashion, resulting in a deflection angle of SOgtep, the orbital period. Port, is 
calculated using Gauss- Chebychev quadrature, and a "representative" diffusion angle 
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during a single orbit, 59 orb, is estimated as 

50 orb = (75 ) 5 9 step 

V -'orb / 

where St is the time step. 
2. The star's velocity vector with respect to the encounter reference frame is calculated 

step ■ 



and a variable L2, which represents the remaining quadratic deflection, is set to 66"^ 



3. The star is tested for entry into the loss-cone, and , if this is case, is removed and 
its mass added to Mbh, whereupon the random- walk is terminated. Otherwise, we 
proceed with the next step. 

4. If L2 < 0, the random-walk is terminated. The star's position and velocity are reset 
to its values before the random-walk procedure. 

5. A new random- walk step is carried out with amplitude A = max [SOorb, min (O.Itt, Agafe, 
and a random direction on the velocity sphere, where Agafe is set to roughly half 
the angular distance to the loss-cone. This way, A becomes progressively smaller 
down to 66orb when approaching the loss-cone in order to keep the risk of missing a 
disruption minimal. According to the new step size, the direction of the star's velocity 
is changed, and L2 is updated: L2 := L2 — A. The random-walk continues at step 3. 

Although many of our results turn out to agree very well with previously published data, 
there are discrepancies when comparing the disruption rates. One possible reason for these 
differences might be related to fa ct that our code uses a shared timestep scheme. In an 



individual timestep scheme, as in 



Freitag fc Bena (1200 ll ). the timestep is some constant 



fraction of the local relaxation time, i.e., dti = f tr{ri), where / is a constant, and the 
subscript i refers to the individual star. In a shared time step scheme the smallest of these 
dti is chosen for all stars. This results in much shorter time steps for stars further out in 
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the c 



uster compared to an individual timestep scheme. As has been noted by 



Freitag et al. 



( 120061 ). the timestep size must be chosen small enough in order to achieve good agreement 



with A^-body simulations, with f < 0.01. While choosing such a small time-step was 



still feasible in the code of 



Freitag fc Benzl (120021 ). to enforce such a criterion for all stars 



in our code would lead to a dramatic slow-down and notable spurious relaxation as the 
timesteps for the stars in th e outer cluster regio ns relative to the local relaxation time 



become extremely small (see 



Giirkan et al. 



2004J ) . In order to reduce the effect of spurious 



relaxation we are forced to choose a larger dt resulting in a larger / for the inner regions, 
up to / ~ 0.1. Fig. [1] shows an example of dt/tr{r) as a function of radius. One can see 
that in this case only for r > O.lr/^ we have dt/tr{r) < 0.01 , while for the inner region it 
rises quickly to 0.1. As a result, any expansion of the inner cluster region, which limits 
the growth of the IMBH in our comparison calculations, will be slower and the disruption 
rates, therefore, larger, due to the higher core densities. Such an expansion occurs either in 
response to the growing IMBH mass or as an initial expansion of the cluster. Previously, we 
applied a procedure that tries to compensate for the larger time steps in the inner region 



(lUmbreit et al. 



20081 ) by evolving the inner stars individually on smaller timesteps while 



keeping the cluster potential consta nt during one shared t ime st ep, a procedure that bears 



some resemblance to the method of 



Marchant &: Shapiro! (1l980l ). However, it turned out 



that in order to achieve good agreement with direct A^-body simulations one has to choose 
the timestep sizes for each cluster configuration separately. This is not only undesirable 
but might also imply that there are additional processes at work that significantly influence 
the disruption rate and are not included in the Monte Carlo scheme. In §4.2.21 we discuss 
several possibilities, among them the wandering of the IMBH. As there is no obvious way 
to compensate for such a processes in a uniform and consistent way through adjustments in 
the two-body relaxation time scale, we do not use the sub-timestep scheme for the present 
paper. 
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w/o BH (Freitag & Benz 2001) 



with BH (Freitag et al. 2006) 
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Fig. 1. — Snapshot of the rat io of global MC 



individual timestep scheme of 



timestep to local relaxation time. While in the 



Freitag fc Bend (120021 ) dt is a constant fraction of tr{r) (dashed 



lines), in our shared timestep scheme dt/tr{r) (solid line) is decreasing with increasing r as 
tr increases. As dt must be chosen large enough to avoid artif i cial re laxation in the shared 



timestep scheme, / is larger for r < than in 
expansion of the cluster. 



Freitag et al. 



(120061 ). resulting in a slower 
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3.2. Initial Conditions 



Table [T] summarizes the initial conditions and main results for all our single-mass 
e [2] fo r our multi-mass runs that include stella r evolution (implemented by 



runs and Tab 


e [21 for our 


Freeeau et al. 


2009 


using 



using the SSE code of 



Hurley et al. 



2002) 



The single-mass clusters consist of stars all of mass 1 Mq with positions and 
velocities chosen according to a Plummer model or a King model with dimensionless 
potential Wq = 10. Radii and times are given in terms of the virial radius, Rvir, and the 
crossing time, tcr, respectively, which are defined by 



Rv 



and 



tr 



5/2 



AEr 



\3/2 ■ 



(4) 



(5) 



where Eq is the total gravitational energy of the cluster. The cluster was evolved up to a 
time Tend, with an initial IMBH mass Mbh,! and rt a constant for all stars. Also shown 
are the resulting final IMBH mass, Mbhj and the number of tidally disrupted stars, Ndisr- 

The initial conditions for multi-mass clusters are chosen as in BMH, with stellar masses 



drawn from a Kroupa mass function (IKroupa 



200l|) in the range of 0.1 - 30 Mq, stars 



evolved at a metallicity Z = 0.001, and initial positions and velocities chosen according to 
a King model with Wq = 7. The tidal disruption radius for each star is given by equation [3] 
with stellar mass and radius provided by the stellar evolution code. In addition to Mbh,!, 
Mbhj and Ndisr the half-mass radius, r^, cluster mass, Mc, and the final relaxation time 
at the half- mass radius, trh, are shown, where the added subscripts i and / indicate initial 
and final values, respectively. 

For each set of initial conditions 9 Monte Carlo runs were performed with varying 
random seed, and the values given in the table are the averages and standard deviations of 



Name 



BME-1 



N 



80, 000 



Wn 



10 



n [Rvir] 



1 X 10- 



MBH,^ [Mc 



0J 



266 



BH,f 



oj 



1429 ± 35 (827) 



END 



[tar] 



3000 



1163 ±35 (561) 



BME-2 



80,000 



10 



1 X 10" 



^00 



1690 ±20 (1388) 



2000 



890 ± 20 (588) 



BME-3 



80,000 



10 



1 X 10- 



2660 



3434 ± 20 (3285) 



2000 



774 ± 20 (625) 



BME-16 



178,800 



10 



1 X 10- 



461 



2171 ±30 (1368) 



2000 



1710 ±30 (907) 



ASFB-50 



100,000 



Plummer 



2.26 X 10- 



50 



10500 ± 100 
(7450, 13050) 



5.4 X 10^ 



10, 450 ±100 
(7.4 X 10^, 1.3 X 10^) 



ASFB-500 



100,000 



Plummer 



2.26 X 10"^ 



500 



11500 ± 100 
(8900, 14500) 



5.4 X 10^ 



11, 000 ± 100 
]A X 10^ 1.4 X 10^) 



able 1 : Details of t he performed Monte Carlo ru ns fo r single-mass clu s ters w ith initial conditions as in 



(j2004a|) (BME) and 



Amaro-Seoane et al. 



f!2004h and 



Baumgardt et al. 



Freitag fc Bend (120021 ) (ASFB). Values in parent 



leses are results 



from the corresponding literat u re. W here two values are given the first one refers to 



second to 



Amaro-Seoane et al. 



fl2004f ). 



Freitag fc Bend (120021 ) and the 



Name 


N 


Wo 


Mbh,, [Mq] 




Mbhj [Mq] 




Mcj [Mq] 


Tend [Gyr] 


^hj [pc] 


\ogtrh [yr] 


BMH-1 


131,072 


7 


125 


4.91 


133 ±3 
(137) 


4± 2 


45,671 ±6 
(45, 534) 


12 


12.02±0.05 
(12.31) 


9.86 
(9.82) 


BMH-2 


131,072 


7 


250 


4.91 


260 ±4 
(280) 


6±2 


45,677±14 
(45,311) 


12 


11.98±0.04 
(12.60) 


9.86 
(9.84) 


BMH-3 


131,072 


7 


500 


4.91 


513 ±5 
(531) 


9±3 


45,400±30 
(44, 741) 


12 


12.36±0.04 
(13.70) 


9.88 
(9.89) 


Table 2: Details of the performed Monte Carlo runs for multi-mass clusters with initial conditions as in 


Baumgardt et al. 





(120051 ) (BMH). The star masses were chosen according to a 



evolution was modeled using the SSE code of 



Hurley et al. 



Krou pa mass function ranging from 0.1 — 30 Mq. Stellar 



mm . 
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these runs. 



4. Idealized Models 

4.1. Imprints of IMBHs 

In Fig. [2] density and velocity profiles from our simulations of single-mass clusters are 
shown, together with the expected r^. calculated from equation ([1]) and using a ^ 0.55 
in A^-body units, appropriate for a Wq = 10 King model. As can be clearly seen, the 
density profile of the inner region of the evolved clusters follow closely the expected 
n{r) oc r~^/^ power-law and the extent of the cusp matches that of the region where the 
velocity dispersion is Keplerian, which in turn matches the expected r^. However, contrary 
to what is seen in direct A^-body simulations, the cusp extends down to much smaller radii, 
especially for black hole masses below 1% of the cluster mass. This is mainly because, in our 
simulations, the central IMBH has a fixed position, while in direct A^-body simulations it is 
allowed to move freely. As a consequence, the density profile flattens inside its wandering 
radius compared to a pure cusp profile, resulting in fewer stars in the central region. As will 
be shown later, this might have an influence on the rate at which stars are tidally disrupted. 



4.2. Disruption Rates 



1 r 



4-2.1. Comparison with \Amaro-Seoane et al.i h2004) and 



Freitag & Beni ^200i 



Fig. B] compares the growth of the IMBH in our simulations and those of 



Amaro-Seoane et al. 



(120041 ). using a gas code, and iFreitag fc Benzl ( 12002| ). using an 
individual time step Monte Carlo code. Here the evolution of a cluster with 10^ stars 



all of 1 Mq and a fixed massive BH in the center with a mass of 50 Vip) or 500 Mp, was 







-20- 




0.01 1 0.01 1 100 

r 

Fig. 2. — Density, p(r), (left panels) and velocity dispersion profiles, cr{r), (right panels) 
for runs BME-1 to BME-3 and BME-16 from top to bottom, respectively, p is given in 
units of Mc/R^j^^ and cr in units of R^ir/tcr- The dotted line indicates the radius of influence 
of the IMBH and the dashed lines represent the theoretically expected power-laws scalings 
p j--'^/^ or (7 ~ r~^/^. 



0.01 1 100 

t/t 

r 



Fig. 3. — BH mass as a function of time for two different initial B H masses, 50 Mrr-, (lower set) 



and 500 Mq (upper set). Shown are results from t he gas code of 



(dash-dotted lines) and the Monte Carlo codes of 
ours (solid lines). 



Amaro-Seoane et al. 



Freitag fc Bena ( 12002| ) (dashed lines) and 



fl2004 ) 
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calculated. The stars were initially distributed according to a Plummer density profile 
with Th = 0.707 DC ■ As can be seen, similar to the results of the Monte Carlo code of 



Freitag &: Bena (120021). our re sults match qualitatively the results of the gaseous method 



of 



Amaro-Seoane et al. 



(120041 ) inasmuch as there is a steep rise of Mbh at the time when 
the cluster densities are largest, which coincides with the time the cluster would formally 
go into core collapse if there were no IMBH. After that time the BH growth levels off due 
to the expansion of the core, described in ^ leading to the convergence of the IMBH 
mass to an asymptotic value. Quantitatively, however, there are differences in the onset 
of the rap id growth phase as well as in the value of the final IMBH mass. The IMBH 



masses m 



Amaro-Seoane et al. 



(120041 ) are generally larger at late times, while the rapid 



growth phase is somewhat de 



ones by 



ayed. Our results follow more closely, and unsurprisingly, the 



Freitag fc Benzl (120021 ) up until shortly after the onset of the rapid growth phase. 



whil e they converge t o larg er Mbh at late times. Part of the reason for this discrepancy 



with 



Freitag fc Benzl (120021 ) must be related to the larger timestep size compared to the 



local relaxation time in the inner region of the cluster and, thus, their slower expansion 
(discussed in §3])- Despite these differences the asymptotic IMBH masses differ by less than 
a factor of two and are, therefore, in reasonable agreement with each other, given the very 
long integration time. 



4-2.2. Comparison with BME 

We now compare the growth rate of IMBHs in our simulations with the direct A^-body 
simulations of BME. For this comparison we restrict ourselves to runs with a larger number 
of particles to ensure that the central cusp around the IMBH is sufficiently populated with 
stars for the Monte Carlo method to be applicable (see Table [1]). 

Fig. m shows the disruption rate as a function of time for runs BME-2 and BME-16 for 
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our Monte Carlo and the direct A^-body code. The quahtative behavior, i.e. the decrease 
of the disruption rate due to the expansion of the cluster is well reproduced. However, our 
rates are systematically larger, always leading to IMBH masses that are larger than in direct 
A^-body simulations (see Tabled]). We find the largest discrepancies for the lowest Mbh/Mc 
ratios, where the total number of disruptions differ by a factor of 2, whereas the disruption 
rates differ by approximately a factor of 1.5. The difference in the disruption rates becomes 
quickly smaller for larger IMBH masses, so that for Mbh/Mc = 1% (Fig. H] (left panel) 
there is agreement within the error bars. However, our rates are still systematically larger 
for that case , leading to approximately 50% more disruptions at the end of the simualtion 
compared to A^-body simulations. This difference is again smaller for Mbh/Mc ~ 3% going 
down to only 25%. 

There are several possible reasons that may explain larger disruption rates in our 
simulations compared to direct A^-body simulations. First, one has to note that our results 
in Fig. m are the averages o yer 9 runs with dif f erent i nitial seeds to generate the cluster, 
while the disruption rates of 



Baumgardt et al. 



(l2004al ) come from only one realization of a 



cluster and are time averages. As the run-to-run variations seem to be larger than the time 
variations in our MC runs, it might be possible that this also applies to the A^-body results, 
in which case the differences in disruption rates could be statistically less significant or even 
disappear for larger BH masses. 

Furthermore, both calculations start from non-equilibrium cluster configurations, 



consisting of a central massive IMBH and a flat constant density core (IBaumgardt et al. 



2004al ). As the Monte Carlo method assumes that the cluster is always in dynamic 
equilibrium, it cannot adequately model the initial phase until such an equilibrium is 
reached. It might, therefore, be possible that the difference in the number of tidal 
disruptions is at least in part due to differences in modeling the initial, violent relaxation 
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process. On the other hand, the number of disruptions during the first 100 tcr account for 
at most 5% of the total in all runs, indicating that these differences have a rather minor 
infiuence on our results. 

Another effect is the wandering of the IMBH due to close passages of stars that are not 
bound to the black hole. Although, one would intuitively think that a wandering IMBH 
would increase the cross section for stars to ent er the tidal radius as it covers a larger 



volume and thus provides a larger c ross section ( 



tends to flatten the density proflle (iBaumgardt et al. 



Chatter iee et al 



2004a 



20021). the wandering also 



Amaro-Seoane et al. 



200J). 



As the disruption rate is proportional to the density near the IMBH, either at Vcru if the 
loss-cone is empty, or if the loss-cone is full, such a flattening of the dens i ty pro flle could. 



Magorrian fc Tremaind ( 1l999l ) discuss 



therefore, cause a lower disruption rate. Similarly, 
the influence of BH wandering on the disruption rate, arguing that, depending on cluster 
parameters and BH mass, the disruption rate can be increased or decreased relative to the 
disruption rate of a flxed central BH. 

In order to estimate this effect more quantitatively we assume for simplicity that, in 
the case of a wandering IMBH, the IMBH remains inside the cusp and the stellar density, 
n{r), outside of the wandering radius, r^, is given by Ucusp ~ r"^/^, while inside it has a 
constant value of ricuspirw)- We furthermore make the assumption that, in this case, we 
are in the full- loss-cone regime, i.e., the loss-cone is constantly replenished with stars on 
a dynamical timescale, contrary to the empty-loss-cone regime, where this happens on a 
relaxation time. The assumption of a full loss-cone for a wandering IMBH is justifled as 
long as its motion is fast enough so that the influx of stars is always sufficiently high to 
completely replenish the loss-cone. The situation is of course different for a fixed BH where 
the loss-cone is not replenished but rather depleted on such a short timescale. In effect, we 
are comparing the disruption rate of fixed BH in the empty loss-cone regime with that of a 
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wandering BH in the full loss-cone regime. 



The disruption rate for a wandering IMBH in the full loss-cone regime is simply given 
by the product of the disruption cross section of the IMBH, T^bh, the stellar number 
density, n*, and the relative velo city dispersion between IM BH and stars, Vrei, evaluated at 



the wandering radius, 7^^, as ^see 



Binney fc Tremaine 



2008 



their eq. 8-123) 



where T^bh is given by 



GM 



BH 



(6) 



(7) 



and Vrei = a/ {v"^ + v'^bh)/'^ with t>* and vbh denoting the stellar and the IMBH velocity 
dispersion, respectively. Using vbh ^ v^, v^{r) = ^jGMBujr in the cusp region, ^ rf,we 
obtain 

n^= ^\fTineVcriy-^\ rt (8) 
where the subscript c refers to the core, Uc = n{ri), and Vc = viji). 

For a fixed IMBH the disruption rate is give n by the influx of s tars into the region 



inside rent-, and can be similarly estimated using fiFrank &: Rees 



mm 



riic 



nircrit) ^BHircrit) vJrcrit, 



(9) 



where = ^ is the loss-cone angle, rc^it the critical radius, and and are replaced 
with Tcrit in equation (JTj). The term 6f^/2 represents the fraction of stars with velocity 
vectors pointing into the loss cone for an isotropic velocity distribution. Using the same 
substitutions as above, equation (Q can be written as 



nic 



5/4 



From the ratio 



nic 



12 



(^crit \ 
Tu, J 



5/4 



(10) 



(11) 
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as well as from equation ([8]) we see, as we expected, that, in general, the disruption rate 
is decreasing with increasing wandering radius. Furthermor e, when we calculate thi s ratio 
using ^ 4 X 10"'^, and Vcrit ~ 2 x 10^'^ from Fig. 1 in (IBaumgardt et al.ll2004al ) for 
a cluster with Mbh ~ 700 Mq and = 8 x 10"^ (run BME-1) we find a value of ^ 0.6 . 
This is very similar to what we obtain for the ratio of disruption rates in BME-16 (see 
Fig. H]). However, equation fllll) does not seem to hold for larger IMBH masses. If we 
calculate the ratio for, e.g., the cluster with Mbh ~ 1395 and = 8 x 10^ in the same 
figure (run BME-2), equation ffTTl) predicts that the disruption rate for a fixed IMBH 
should be lower than for a wandering IMBH by a factor of ~ 1.6 while from Fig. Hlwe find 
that it is mostly larger. The reason for that is probably related to the assumption of a full 
loss-cone when deriving equation ([8]), as it is clear that in the limit of very small r^, and 
consequently larger IMBH masses, the empty loss-cone regime, and thus, equation (ITUI) . 
must be approached. Indeed, in Fig. H] (left panel) we see that the disruption rates seem 
to converge at late times. Similarly, we find very good agreement in the total number of 
disruptions between A^ -body and our Monte Carlo runs for IMBHs with even larger masses. 



As r„ 



M^/Mbh (jChatterjee et al. 



20021 ) it, therefore, appears that the effect of a 



wandering IMBH has a negligible influence on the disruption rate for Mbh/M^, > 1000 

A third effect that could decrease the disruption rate and that is not included in the 
current Monte Carlo scheme is the formation of, and the strong interaction with, an IMBH 
binary that is able to scatt er other stars into the outer cluster regions. This mechanism has 
been recently sugge sted by 



Gil 



cluster with IMBH. iGill et al. 



et al. 



(120081 ) to suppress mass-segregation in a multi-mass 



(120081 ) also found that the efficiency of this scattering process 



is relatively independent of the relative mass of the IMBH. Given the trend we see for the 
difference in the number of tidal disruptions between our simulations and A^-body results as 
a function of Mbh/Mc, it appears that this mechanism is unlikely to be the main reason. 
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Finally there is the possibihty that, as aheady argued in §4.2.11 the larger disruption 
rates are a result of the slower expansion of the inner regions caused by the larger timestep 
size relative to the local relaxation time. However, while this slower expansion might 
increase the disruption rate in general, it is not at all clear why this should make a 
larger difference for clusters with lower-mass IMBHs. For instance, when we compare the 



ratios of the fina^ 



masses from our runs to the results of the individual time step code of 



Freitag &: Bena (|2002| ) for the different Msu^i we see from Fig. [3] that they are rather 
similar (0.70 for Mbh^i = 50 Mq and 0.71 for MBH,i = 500 M©). On the other hand, the 
final mass ratios for runs BME-1 to BME-3 vary dramatically in comparison (from 0.6 to 
0.95) for a similar range in initial IMBH masses. Another reason why it is less likely that 
the differences to the direct A^-body results are caused by differences in the expansion rates 
is given in Fig. where we compare the Lagrange radii for run BME-16. As can be seen, 
apart from very minor deviations, the expansion of the cluster in our simulation agrees well 
with the A^-body results, even for the regions within O.lrh- However, one should also note 
that, here, only the radii are shown that contain more than 1% of the total cluster mass, as 
only for those A^-body data were available. Larger differences would be expected for smaller 
radii, especially because for this particular run the 1% Lagrange radius is just outside the 
cusp region which mostly determines the disruption rate. 

In general we can say that, despite the significant differences in the total number of 
disrupted stars, the disruption rates do not differ greatly from the ones from direct A^-body 
simulations and are even in very good agreement for runs with Mbh/Mc > 0.01 given the 
error bars. Furthermore, from Fig. [5] we see that our Monte Carlo code can reproduce the 
evolution of the cluster structure, at least down to the cusp region, quite well. We expect 
that this also applies to other clusters as long as the wandering radius of the IMBH is inside 
the cusp. 
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5. Realistic Cluster Models 

5.1. Comparison to BMH 

Here we re-examine a subset of the A^-body models of BMH using our Monte Carlo 
code to see if we are able to reproduce their results and, in particular, the inner surface 
brightness slopes mentioned in §2J For this comparison we only include runs with large 
in order to populate the high mass end of the IMF sufficiently well, which is important for 
the applicability of the Monte Carlo method. Furthermore, as in BMH, we only include 
bright stars, defined as main-sequence stars and giants with masses larger than 90% of the 
turn-off mass, for the calculation of the surface brightness profiles to take into account that 
those stars contribute most (> 80%) of the observed light. 

Table [2] gives a summary of the resulting cluster structural parameter after the clusters 
have been evolved to an age of 12 Gyr together with the corresponding results of BMH. In 
contrast to our single-mass runs we find lower IMBH masses in our runs compared to the 
direct A^-body results. However, this difference is much less pronounced, if not negligible, 
as it is at most 20 Mq, resulting in final IMBH masses that deviate at most by 7% from 
the A^-body results. The reason for that is not only because the tidal disruption radii of 
the stars are generally lower compared to our single- mass cases but, more importantly, the 
cluster had a much lower density in the core initially and subsequently expanded because 
of mass loss due to stellar evolution. After the core started contracting again, after about 
1 Gyr, it became quickly dominated by massive dark remnants, mostly black holes and 
massive white dwarfs, that drove out lower-mass main-sequence stars. This decreases the 
overall disruption rate, because the main-sequence stars which are more easily disrupted 
given their much larger radii compared to compact remnants, are, on average, much further 
out, while the compact remnants, though much closer to the IMBH, are unlikely to get 
disrupted. 
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The fact that for these simulations there are fewer disruptions in our Monte Carlo runs 
than in the A^-body simulations is most probably again related to the wandering of the 
IMBH. While for a fixed IMBH there is a well populated cusp, for a wandering IMBH no 
such clear cusp can be identified. In the A^-body simulations the massive dark remnants in 
the cusp are ejected through strong gravitational interactions with the IMBH, which allows 
the main-sequence stars and giants to diffuse inwards and to come closer to the tidal radius. 
In addition, due to its motion, the IMBH is also able to come closer to stars that are just 
outside of the cusp region, which also increases the number of disruptions. 

Comparing the cluster structural parameters we find again very good agreement 
between the Monte Carlo and direct A^-body runs. There are only minor but systematic 
differences inasmuch as in our Monte Carlo runs the clusters have systematically larger 
masses, though only by less than 1.5%, and are more compact, with their half-mass radii 
differing at most by 10% . The most likely reason isthat, in our simulations, we do not 
model close encounters with stars tightly bound to the IMBH. Those interactions might 
frequently lead to the ejection of massive stars or remnants, which also contributes to the 
cluster expansion. In addition, the somewhat larger disruption rates might have also added 
to the stronger expansion. However, given the minor differences between the Monte Carlo 
and A^-body results, close encounters with stars tightly bound to the IMBH do not seem to 
have a significant infiuence on the cluster structure as a whole. 

We also achieve good agreement for the surface density profiles between the two 
methods. FigJHl shows the two-dimensional density profiles of bright stars for two clusters at 
the end of the simulation. These profiles are obtained by averaging five snapshots obtained 
at 50 Myr intervals. The profiles show only very shallow cusps with a power-law slopes a 
between —0.2 and —0.3, consistent with the A^-body results. In order to reduce the amount 
of noise and obtain a more reliable fit of the inner slope of the surface brightness profile, we 
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re-ran run BMH-2 but with twice as many particles. The resulting profile is shown in Fig. 



El As one can see, t 
a = —0. 25 found in 
found by Miocchil feoOTf ) . 



le power-law fit has a s lope of a = —0.23, very close to the average 



Baumgardt et al 



(120051 ) and agreeing very well with the relation 



5.2. Comparison to Real Star Clusters 



The ultimate goal of any cluster simulation is to reproduce observations of real star 
clusters as closely as possible. For globular clusters these observations are mostly in the 
form of photometric data and surface brightness profiles, while there are only a few clusters 
for which well measured velocity dispersion profiles are available. In the previous section 
we found that our Mo nte Carlo simulati o ns wi th IMBH as well as the corresponding 



A^-body simulations of 



Baumgardt et al. 



(120051 ) are able to reproduce the intermediate 



inner surface brightness slopes seen in some Galactic globular clusters. As discussed in 
^ an IMBH also influences the surface brightness profile in that it produces rather large 
cluster cores as measured by Vc/rh- It is, therefore, interesting to see if the combination 
of slope and concentration we find in our models matches any observations. However, 
instead of comparing slopes and concentrations quoted in the literature it is more suitable 
to directly compare surface brightness profiles in this case, given that the tidal ra dii of 



observed clusters are very uncertain (see, e.g., discussion in iBaumgardt et al. 



20091). 



For this comparison we choose the cluster NGC 5694 for which 



Noyola fc Gebhardt 



(120061 ) report an inner slope of 0.19 ± 0.11, close to our fit in Fig. [6] (right panel). In 
addition, this cluster is particularly well suited as it is sufficiently isolated, with a distance 
of ~ 30kpc from the Galactic center. This is because all our comparison calculations, 
as well as the corresponding A^-body simulations, do not include tidal effects. Since our 
profiles represent projected luminosity profiles (approximated by the number of bright stars 
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per unit area), while the observed ones are given in magnitudes per square arc seconds, we 
first converted magnitudes to luminosities. We then applied linear scaling factors to the 
radius and brightness in order to match the observed profile to our our results as closely as 
possible. 

Fig. [7] shows the resulting profile together with our profile from Fig. O As one can see, 
we obtain an extremely good match within 2rh, while outside, the observed profile seems 
to fiatten. This excess light is usually attributed to both the confusion with background 
stars and stars that are esc aping but remain still close to the cluster for a considerable time 



( iFukushige fc Heggie 



2OOOI ). which we cannot model using a simple tidal cut-off prescription. 
In this case, however, confusion with back ground stars is a more likely explanation as the 



estimated tidal radius from the literature ( iHarris 



19961) is larger by a factor of 4 compared 



to the radius where the observed profile begins to fiatten. In addition, when we correct 
for the background by simply subtracting the value of the outermost point from the other 
points in the profile, we reach an even better agreement with our profile out to 3 r^- 



Miocchil (120071 ) ■ NGC 



This close agreement is somewhat surprising as, according to 
5694 is not expected to harbor an IMBH based on the fact that its concentration is too 
large for the relatively steep inner surface brightness slope (see relation [2]) . Howe ver, 



this s tatement is based on values for the concentration quoted from the literature ( IHarris 



19961 ) which, as we mentioned before, are in general rather uncertain due to difficulties 
determining the tidal radius of a cluster. An error of 40% in the measured tidal radius 
would bring the value of the slope and the concentration of the cluster in agreement with 
relation ([2]). Furthermore, the observed slope has a rather large uncertainty as well, which, 
even considered alone, could make the observed parameters consistent with this relation. 
Therefore, it appears that the close match we obtain for the shape of the surface brightness 
profile in our simulation with the observed one of NGC 5694 does not necessarily contradict 
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the validity of the slope-concentration relation found by 
broadly consistent with it. 



Miocchil (120071 ) but seems to be 



6. Summary and Conclusions 

In this paper we studied the influence of IMBHs on the evolution of globular clusters 
using our Monte Carlo code, which has been extended to include the effects of a central 
IMBH on the stellar distribution. The IMBH is treated as a fixed point mass in the cluster 
center, and at each time step stars are tested for entry into their loss-cone. In order to test 
our implementation we carried out a large number of idealized, as well as more realistic 
dynamical simulations of globular clusters and compared our results with published results 
from direct A^-body and Fokker-Planck calculations. 

In general we found that our results agree reasonably well with results from A^-body 
and Fokker-Planck type codes. Significant differences only exist for idealized models while 
for more realistic clusters these differences become negligible. 

We found that our code is able to reproduce the expected density cusp and Keplerian 
velocity profile around the IMBH for a single-mass cluster very well. In addition, we also 
found that the evolution of the cluster density profile outside of the 1% Lagrange radius is 
in very good agreement with direct A^-body simulation. The only notable difference is that 
the cusp extends f urther down towards the IMBH than in the direct A^-body simulations of 



Baumgardt et al. 



(j2004al ). The reason is most likely that the IMBH is fixed at the cluster 
center while in the A^-body simulations it is allowed to move freely. Although this has 
a minor effect on the density profile, causing it mainly to be flatter inside r^, we argue 
that it produces significant differences in the disruption rates for single-mass clusters. In 
particular, as we demonstrated by simple estimates of the disruption rate, the effect of 
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IMBH wandering might be able to explain why we see up to a factor of 2 larger disruption 
rates for low-mass IMBHs than for higher-mass ones, compared to A^-body results, which 
is difficult to explain otherwise. For a sufficiently fast moving IMBH the loss-cone can 
always be assumed full in which case the disruption rate is a strong function of the density 
at the disruption radius rj. Therefore, for larger r^, and consequently for lower Mbh, the 
density at becomes lower inside the very steep power-law cusp which, in turn, decreases 
the disruption rate. In our calculations the disruption rate is determined by the diffusion of 
stars into the loss-cone. This rate turns out to be larger than for the case of a wandering 
IMBH for our cluster parameters if is significantly larger than Vcrit, which agrees with 
the differences we see for run BME-16. For lower r^j, and, thus, larger Mbh, loss-cone 
depletion effects appear to become important as we find from our simulations that the rates 
for these two cases converge for IMBHs with Mbh/M^ > 1000, while equation (fTTi) predicts 
that the rate for a wandering IMBH should become larger than for a fixed one. 

Comparing the IMBH mass growth in the em pty-loss-cone regime with d ifferent 



Fokker-Planck type c odes, such as th e 



code by 



Amaro-Seoane et al. 



(120041 ) and the 



Monte Carlo code by iFreitag fc Benzl (j2002l ). we found better agreement with our results. 
For all codes the IMBH growth shows the same qualitative behavior, that is, the the 
masses converge to a constant value. Only at late times the masses differ significantly by 
factors of 1.3 to 2. Th e larger IMBH masses in our simulations compared to the results of 



Freitag &: Bend (120021 ) are most likely due to the slower expansion of the inner regions in 



our simulations, which are caused by the larger timesteps relative to the local relaxation 
time, a consequence of the shared timestep scheme used in our code. Contrary to the effect 
of IMBH wandering, this effect does not introduce a dependence on Mbh and our final 
masses always differ by 30% regardless of the initial Mbh- 



The situation for realistic cluster models is, however, rather different. In this case the 
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disruption rates are generally lower because the cluster core gets quickly dominated by dark 
remnants, which have extremely small r^. Therefore, the influence of the IMBH wandering 
is, in general, much reduced and does not produce signiflcantly different results compared 
to a flxed IMBH. There are only minor deviations in the flnal cluster mass and size, with 
our clusters being more compact by ^ 10% and slightly more massive by ~ 1%. The reason 
here might be related to the fact that we do not model close encounte rs with stars in th e 



cusp, which could be an important ejection mechanism for the stars (iBaumgardt et al. 
2004al ). S uch encounters are also responsible to efficiently remove dark remnants from 



the cusp ( iBaumgardt et al. 



2004bl ) and, thus, allow for lower-mass main-sequence stars to 
get close to the IMBH and disrupted. The somewhat larger number of disruptions in the 
A^-body simulations is also caused by the wandering of the IMBH as it can, this way, get 
closer to the main-sequence stars that are, due to mass segregation, further out. Although, 
the total mass the IMBHs gain in the A^-body simulations is up to 3 times larger than in 
our simulations, the total difference amounts to just 20 M0 at most or less than 7% of the 
total cluster mass and is, thus, very minor, given that this difference comes from not much 
more than 10 disrupted main-sequence stars. 

Using our Monte Carlo code we are also able to reproduce the intermediate power-law 



slopes in the sur: 



Baumgardt et al 



ace br ightness profile seen in all A^-body simulations with IMBH in 



(I2OO5I ). Although, there is considerable scatter in the profile which 
makes the determination of the individual slopes rather uncertain, we show that they, 
nevertheless, are within the same region as in the A^-body simulatio ns. Using twice as many 



stars in our simulations than were used in 



Baumgardt et al. 



(120051 ) we have a high enough 



resolution to obta i n a re liable fit, which turns out to be close to the average slope found by 



Baumgardt et al 



(I2OO5I ). This, once again, confirms that, despite differences in the details 
of the dynamics in the innermost cusp region, the overall cluster structure and evolution is 
rather well reproduced. 
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In order to see if our cluster models can reproduce observations, we calculated the 
surface brightnes s profile from our largest N run and compared it to the profile of NGC 



5694 obtained by 



Noyola &: Gebhardt 



(I2OO6I ). a cluster similarly isolated from the Galactic 
tidal field as our models. Here we found that the overall shape of the cluster matches the 
observed profile up to 2 half-mass radii very well. Outside of that region the observed profile 
becomes flatter which is most likely due to confusion with the background. Although, the 
existenc e of an IMBH in this particular cluster seems inconsistent with the static models of 



Miocchil ( 120071 ). with NGC5694 having a concentration that is too large for the measured 
inner surface brightness slope, considering; possible errors in the value for the concentration 



quoted in the Harris catalog ([Harris 



19961 ) and the rather large uncertainty in the inner 
surface brightness slope, this is not necessarily the case. The fact that our model with 
IMBH not only reproduces the inner slope but also matches the shape of the observed 
surface brightness profile of NGC 5694, up to a radius of 2 r/j, provides a stronger indication 
for the existence of an IMBH with a mass of the order of 0.3% of the total cluster mass. 

In summary, we find that our Monte Carlo code, which now includes the effects of a 
central IMBH on the stellar distribution, is able to reproduce the overall structure and 
evolution of single-mass as well as of more realistic, multi-mass cluster models reasonably 
well compared to direct A^-body simulations. There are differences in the disruption rates 
and IMBH masses which are most likely related to the wandering of the IMBH, which is 
not included in the Monte Carlo code. However, these differences become almost negligible 
for realistic cluster models as the disruption rates are generally very low for these cases. 
In addition, since the 



(IChatterjee et al. 



MBH wandering radius scales with Mbh/M^ as ~ ■>/ Mbh/M^ 
2OO2I ) the influence of the IMBH motion is reduced for more massive 
clusters with a similar Mbh/Mc- One process that still plays a role for more massive 
clusters and which is not included in our Monte Carlo simulations are the close gravitational 
interactions with cusp stars which we plan to implement in the near future. 
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Fi g. 4. — Comparison of the disruption rates per crossing time for the direct A^-body results 



of 



Baumgardt et al. 



(j2004al ) (open symbols) and for our MC simulations (filled symbols). 
Shown are results of the runs BME-2 (left) and BME-16 (right). The disruption rates of 
the MC runs are the average rates from 9 runs with different random seeds, and the error 
bars are the standard deviation of this average. The disruption rates and error bars of the 
A^-body results are time averages and their standard deviations respectively. 



- 38 - 



> 



1 1 1 — I — I — I I I I 1 1 1 — I — I — I I I I 1 1 1 — I — I — I I I I r 

- - - - - , ;^ , y ,, r^ , ni ,. TnT i W l W i ffT iiWl l l l l lB i lM l ff l» »i li ag 



MuiiigiiriijBiTiiniiiiifrgfl' 



- Tr , B . -7 i|i riHiWF >i ?Wii;i li H )| f ili 



- rif ' r i ri ). T»i<iffi i ffi ii a > w »l|iiO 



I 




Fig. 5. — Evolution of the Lagrange radii of run BME-16 for our Monte Carlo simulations 
(solid lines) and direct A^-body simulations (dashed lines). Shown are the radii that contain 
1%, 2%, 5%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, and 90% of the total mass of 
the cluster. The cluster expands due to the energy generated in the cluster center by the 
tidal disruption of stars. Apart from a small jump in the A^-body data for the outer 90% of 
the cluste r, which appears t o be sp urious, our results are in reasonable agreement with the 



results of 



Baumgardt et al. 



(!2004ah . 
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Fig. 6. — Surface density profile of bright stars for two clusters with different numbers of 
stars and BH masses. The dashed line in the right panel is a power-law fit to the inner 
region of the cluster, while the two dashed lines in the left panel are power-laws with slopes 
bracketing the range [—0.2, —0.3] suggested by BMH for clusters harboring IMBHs. The 
inner parts of our surface density profiles are in good agreement with the results of BMH. 
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Fig. 7. — Surface brightness profile (SBP) from one of our MC models (thin dashed line) 



and for NGC 5694 fNoyo 



IMBH f lBaumgardt et al 



a fc Gebhardt 2006; filled circles), a cluster that might harbor an 



20051 ). showing a shallow cusp. Also shown is a power-law fit to our 



model in the inner region with a slope of -0.23, very close to the value of -0.19 determined 
from the SBP of NGC 5694 (Noyola & Gebhardt 2006). 
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